#include "include/MOL.h"
#include <ctime>

using namespace std;

// r = k * nu / h^2
double nu = 1;
double r = 1;

// 空间和时间步长
double h = 1.0 / 64;
double k = r * h * h / nu;

// 空间和时间规模
int SN = 64;
int TN = 10;
double T = TN * k;

// 边界条件
Vector g(Real t)
{
    return Vector(SN + 1, 0);
}

// 初值条件
Real u(Real x)
{
    Real r = 0;
    if (9.0 / 20 <= x && x < 0.5)
    {
        r = 20 * (x - 9.0 / 20);
    }
    if (0.5 <= x && x < 11.0 / 20)
    {
        r = -20 * (x - 11.0 / 20);
    }
    return r;
}

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-27 22:10
 *
 * @description: ESDIRK 方法
 */

void testMOL()
{
    // 计算区域
    Domain<1> domain({0}, 1);

    // 时间网格
    TimeGrid G(TN + 1, Vector(SN + 1, 0));

    // 循环填充初值条件
    do
    {
        Real x = (Real)domain.gridPoint(G[0]);
        G[0].value() = u(x);
    } while (G[0]++);

    // 系数矩阵 A
    Matrix A(SN + 1, 0);
    int M = SN + 1;
    for (int i = 1; i < M - 1; i++)
    {
        A[i * M + i] = -2 * nu / h / h;
        A[i * M + i + 1] = nu / h / h;
        A[i * M + i - 1] = nu / h / h;
    }
    A[0] = A[M * M - 1] = -2 * nu / h / h;
    A[1] = A[M * M - 2] = nu / h / h;

    // 定义求解方法
    ESDIRK method(k);

    // 线方法求解器
    MOL mol;
    mol.solve(A, G, g, &method);

    cout << "U = [" << G.back() << "];" << endl;
    cout << "plot(U);" << endl;
}

int main(int argc, char **argv)
{
    testMOL();

    return 0;
}